Für das Supply Chain Management der Daimler AG existiert ein Track und Trace System, das bestimmte Daten aufzeichnet. Ziel ist es aus den Daten nützliches Wissen für die Fachbereiche zu generieren, sodass die Lieferungen besser gesteuert werden können. Im ersten Schritt werden die Daten einer explorativen Datenanalyse unterzogen, um erste Erkenntnisse über die Lieferungen zu bekommen. Im Anschluss werden auf Top-Level Ebene einige Modelle entwickelt, welche die Ankuftszeiten vorhersagen. Dabei handelt es sich um erste Modelle, um die Einsatzmöglichkeiten zu evaluieren. Langfristig werden weiterhin diese Ziele verfolgt: -Forecast der Transportzeiten zwischen den einzelnen Stationen -Demand in den Produktionsstätten (entweder forecasten oder falls vorhanden mit aufnehmen) und anschließend mit aktuellen Lieferungen abgleichen, ob diese rechtzeitig da sind. -Klassifikation von Lieferungen in Risikogruppen
library(tidyverse)
library(ggmap)
library(plotly)
library(dplyr)
library(lubridate)
library(ggridges)
library(ggplot2)
library(RColorBrewer)
library(bupaR)
library(scales)
library(mltools)
library(caret)
library(data.table)
library(mlr)
library(xgboost)
library(Metrics)
library(MLmetrics)
library(keras)
library(randomForest)
library(e1071)
raw_data <- read.csv("~/R 3.5 Files/Fallstudie Logistik/data_sample_case_study.csv", sep = ";")
raw_event_location <- read.csv("~/R 3.5 Files/Fallstudie Logistik/event_locations_case_study.csv", sep = ";")
#Variablen
colnames(raw_data)
## [1] "ContainerNumber" "ShipmentNumber" "PoCreationDate"
## [4] "Material" "ShipmentQuantity" "Controller"
## [7] "EventName" "EventTime" "EventMessageTime"
## [10] "EventLocation" "VesselName" "PortofDischarge"
## [13] "Time2Arrival"
#Überblick der Daten
head(raw_data)
## ContainerNumber ShipmentNumber PoCreationDate Material
## 1 HLXU8170197 4868820 19.02.2018 00:00 A2760101414
## 2 HLXU8170197 4868820 19.02.2018 00:00 A2760101414
## 3 HLXU8170197 4868820 19.02.2018 00:00 A2760101414
## 4 HLXU8170197 4868820 19.02.2018 00:00 A2760101414
## 5 HLXU8170197 4868820 19.02.2018 00:00 A2760101414
## 6 HLXU8170197 4868820 19.02.2018 00:00 A2760101414
## ShipmentQuantity Controller EventName EventTime
## 1 72 T15 Loading Ship 07.04.2018 13:58
## 2 72 T15 Departure Ship 07.04.2018 19:05
## 3 72 T15 Unloading Ship 21.04.2018 06:19
## 4 72 T15 Departure Truck 21.04.2018 11:42
## 5 72 T15 Arrival Truck 08.05.2018 18:12
## 6 72 T15 Arrival Truck 08.05.2018 18:12
## EventMessageTime EventLocation VesselName PortofDischarge
## 1 08.04.2018 03:19 urn:jaif:id:loc:26LNLRTM NYK ROMULUS
## 2 09.04.2018 16:29 urn:jaif:id:loc:26LNLRTM NYK ROMULUS
## 3 21.04.2018 11:59 urn:jaif:id:loc:26LUSCHS NYK ROMULUS
## 4 21.04.2018 13:29 urn:jaif:id:loc:26LUSCHS NYK ROMULUS
## 5 23.05.2018 22:04 urn:jaif:id:loc:26LUSCHS NYK ROMULUS
## 6 23.05.2018 17:32 urn:jaif:id:loc:26LUSCHS NYK ROMULUS
## Time2Arrival
## 1 125
## 2 125
## 3 111
## 4 111
## 5 94
## 6 94
#Summary Shipment Quantity
summary(raw_data$ShipmentQuantity)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.00 24.00 60.00 60.62 84.00 162.00
#Shipment Quantity Boxplot
plot_ly(y = raw_data$ShipmentQuantity, type = "box") %>%
layout(title = "Boxplot Shipment Quantity",
xaxis = list(title = ""),
yaxis = list(title = ""))
Per Geocoding können die Koordinaten der Event Locations gezogen werden. Diese werden für eine spätere Visualisierung notwendig
register_google(key = "AIzaSyCwltxXK3VXCx0SDJxBPdfLMQWcEUOEvGk")
location <- as.character(raw_event_location$Bezeichnung)
cords <- geocode(location)
raw_event_location$lat <- cords$lat
raw_event_location$long <- cords$lon
#Left join, um die Bezeichnung sowie die Koordinaten in den gesamten Dataframe zu bekommen
data <- left_join(raw_data, raw_event_location, "EventLocation")
#Neue Spalte, die anzeigt wenn ein Teil angekommen ist
data$ShipmentArrived <- if_else(data$EventName == "Goods Receipt Dock" & data$Bezeichnung == "Werk Tuscaloosa", TRUE, FALSE)
#Umwandeln der Spaltentypen in Standardformat Zeit
data$EventTime <- as.POSIXct(as.character(data$EventTime),format = "%d.%m.%Y %H:%M")
data$EventMessageTime <- as.POSIXct(as.character(data$EventMessageTime),format = "%d.%m.%Y %H:%M")
data$PoCreationDate <- as.POSIXct(as.character(data$PoCreationDate),format = "%d.%m.%Y %H:%M")
data$PoCreationDate <- as.Date(data$PoCreationDate)
#Als Numeric speichern
data$ShipmentNumber <- as.numeric(data$ShipmentNumber)
#Monat, Jahr, Wochentag und Quartal hinzufügen
data$monat <- month(data$PoCreationDate)
data$jahr <- year(data$PoCreationDate)
data$wochentag <- weekdays(data$PoCreationDate)
data$quater <- lubridate::quarter(data$PoCreationDate, with_year = TRUE, fiscal_start = 1)
#Daten 18.02.2018 bis 23.04.2019
summary(data$PoCreationDate)
## Min. 1st Qu. Median Mean 3rd Qu.
## "2018-02-18" "2018-06-25" "2018-09-25" "2018-09-24" "2019-01-08"
## Max.
## "2019-04-23"
#Daten 05.04.2018 bis 24.05.2019
summary(data$EventTime)
## Min. 1st Qu. Median
## "2018-04-05 14:01:00" "2018-08-19 11:51:45" "2018-11-20 13:16:00"
## Mean 3rd Qu. Max.
## "2018-11-17 20:47:29" "2019-03-03 16:09:15" "2019-05-24 22:10:00"
## NA's
## "4"
#4 Einträge von EventTime haben ein NA
#Vorgehen: Da EventTime und EventTimeSend meistens ziemlich gleich sind werden die NA's damit ersetzt
data$EventTime <- if_else(is.na(data$EventTime) == TRUE, data$EventMessageTime, data$EventTime)
#Testen, welche Spalten weitere NA's haben
colnames(data)[colSums(is.na(data)) > 0]
## [1] "Bezeichnung" "lat" "long"
#Filtern auf NA's
bez <- data %>%
filter(is.na(Bezeichnung) == TRUE)
#Es ist keine Event Location vorhanden -> deswegen gibt es kein Mapping
#Für die genauere Analyse wird eine Sendungsnummer der betroffenen Events exploriert
bsp <- data %>%
filter(ShipmentNumber == 5120761)
#Es sieht so aus als wäre diese Zeile ein Duplikat in dem die Event Location nicht richtig aufgezeichnet wurde
#Aus diesem Grund werden nun alle Duplikate entfernt, basierend auf der ContainerNummer, der ShipmentNumber, des ProCreationDate, des EventName, des Material, des Controller sowie der EventTime
data <- data[!duplicated(data[1:8]),]
#Testen, ob noch NA drin sind
colnames(data)[colSums(is.na(data)) > 0]
## [1] "Bezeichnung" "lat" "long"
#Es sind noch 53 NA's enthalten -> wieder überprüfen
#Filtern auf NA's
bez <- data %>%
filter(is.na(Bezeichnung) == TRUE)
bsp <- data %>%
filter(ShipmentNumber == 5281788)
#Das sind alles Event mit dem Name ContainerClosed -> Das sind alle, die auch ein NA haben. Dementsprechend kann keine eindeutige Zuordnung des Events stattfinden.
#Dies sollte mit dem zuständigen Fachbereich geklärt werden!
#Für jetzt werden diese Events als eigenstädnige Event Location behandelt. Es wird sicher später zeigen, wo die Warenströme nach diesem Event hingehen
data$Bezeichnung <- if_else(data$EventName == "Container Closed", "Container Dummy", as.character(data$Bezeichnung))
#Die restlichen NA's in Latitude und Longitude können erstmal so gelassen werden, da keine Adresse vorhanden ist
#Diese dataframes wieder löschen, um RAM zu sparen
remove(bez)
remove(bsp)
Dieser Plot zeigt die monatlichen Mengen an, die transportiert werden getrennt nach Jahren. Es ist zu erkennen, dass 2019 bisher weniger transportiert wurde als im Vorjahr. Das Maximum lag im Jahr 2018 im April.
#Die Daten auf Basis der Sendungen
df_sendungen <- data[!duplicated(data[2]),]
#Nur auf monatlicher Ebene
df_sendungen <- df_sendungen %>%
mutate(Date = format(as.Date(df_sendungen$PoCreationDate), "%Y-%m"))
#Anzahl an Jahre -> durch dieses Vorgehen wird sichergestellt, dass dieser Plot für eine beliebige Anzahl an Jahren richtig angezeigt wird
Years = as.data.frame(unique(year(df_sendungen$PoCreationDate)))
colnames(Years) <- c("Jahre")
#Plot initialisieren
p <- plot_ly(type = 'scatter', mode = 'lines', line = list(shape = "spline")) %>%
layout(title = "Monatliche Quantity in Einheiten",
xaxis = list(title = ""),
yaxis = list(title = ""))
#Eigenen Farbvektor erstellen
color_vec <- c('rgb(8,48,107)','rgb(205, 12, 24)','rgb(0, 0, 0)','rgb(128, 128, 128)','rgb(0, 255, 255)', 'rgb(218, 165, 32)','rgb(34, 139, 34)')
#Für jedes Jahr wird ein Trace erstellt
i = 1
while (i <= nrow(Years)) {
#Auf Jahr i filtern
df_neu <- df_sendungen %>%
filter(jahr == Years$Jahr[i])
#Quantity auf monatlicher Ebene aggregieren
eda <- aggregate(df_neu$ShipmentQuantity, by=list(ds=df_neu$Date), FUN=sum)
eda$ds <- paste(eda$ds, "01", sep = "-")
eda$x <- as.numeric(eda$x)
eda$ds <- as.Date(eda$ds)
eda$month <- month(eda$ds)
name <- unique(year(eda$ds))
#zum Plot hinzufügen
p <- p %>%
add_trace(data = eda, x = ~month, y = ~x , name = paste("Quantity ",name), line = list(color = color_vec[i]))
i = i+1
}
p
Dieser Plot zeigt die Dauer von jedem Hafen bis nach Tuscaloosa. Hier ist wohl angemerkt, dass die Daten nicht auf der Eventebene betrachtet werden. Das bedeutet, dass unterschiedliche Events die Daten verzerren
#Die gleichen Werte bei unterschiedlichen Bestellungen entfernt
df_p1 <- data[!duplicated(data[4:11]),]
df_p1 <- df_p1 %>%
select(ShipmentNumber, Bezeichnung, Time2Arrival)
#Hier wieder die Duplikate entfernen -> Sonst verzerrt
df_p1 <- df_p1[!duplicated(df_p1[1:3]),]
plot_ly(df_p1, y = ~Time2Arrival, color = ~Bezeichnung, type = "box")
Dieser Plot zeigt die Dauer von jedem Hafen bis nach Tuscaloosa. Hier ist die Detailebene genauer, da die Darstellung auf der Hafen und Eventebene liegt. Es ist deutlich zu erkennen, dass bestimmte Häfen sehr hohe Ausreißer haben (z.B. Bremerhafen oder Charleston)
#Die gleichen Werte bei unterschiedlichen Bestellungen entfernt
df_p1 <- data[!duplicated(data[4:11]),]
#Hafen und Events verbinden
df_p1$HafenEvent <- paste(df_p1$Bezeichnung, df_p1$EventName, sep = "-")
df_p1 <- df_p1 %>%
select(ShipmentNumber, HafenEvent, Time2Arrival)
#Hier wieder die Duplikate entfernen -> Sonst verzerrt
df_p1 <- df_p1[!duplicated(df_p1[1:3]),]
plot_ly(df_p1, y = ~Time2Arrival, color = ~HafenEvent, type = "box") %>%
layout(showlegend = FALSE)
Dieser Plot zeigt die Häufigkeit der verwendeten Schiffe. Die ersten drei Schiffe werden ungefähr gleich häufig verwendet. Das spätere Modell wird zeigen, welchen Einfluss die Verwendung der Schiffe auf die Ankunftszeit hat.
#Häufigkeit der Schiffe
df_p2 <- count(df_sendungen, vars = VesselName)
#Absteigend ordnen und als Factor speichern
df_p2 <- df_p2[order(df_p2$n),]
df_p2$vars <- factor(df_p2$vars, levels = c(as.character(df_p2$vars)))
#Plottet die Häufigkeit der Schiffe
plot_ly(df_p2, x = ~n, y = ~vars, type = 'bar',
marker = list(color = 'rgb(158,202,225)',
line = list(color = 'rgb(8,48,107)', width = 1.5))) %>%
layout(title = "Häufigkeit der Schiffe",
xaxis = list(title = ""),
yaxis = list(title = ""))
Dieser Plot zeigt den Materialbedarf in den einzelnen Quartalen an. Es ist zu erkennen, dass der Materialbedarf im 1. Quartal 2019 leicht über dem des Vorjahres liegt. Mit mehr Daten könnten hier Zusammenhänge besser interpretiert werden.
#Daten nach Quartal gruppieren und die Mengen summieren
df_p3 <- df_sendungen %>%
group_by(quater) %>%
summarise(Quantity = sum(ShipmentQuantity))
#Achsen richtig benennen und den Durchschnitt als Vergleichsgrundlage berechnen
df_p3$QuaterName <- paste("Quartal", df_p3$quater, sep = " ")
df_p3$avg <- mean(df_p3$Quantity)
p <- plot_ly(df_p3) %>%
add_trace(x = ~QuaterName, y = ~Quantity, type = 'bar', name = 'Materialbedarf',
marker = list(color = 'rgb(158,202,225)',
line = list(color = 'rgb(8,48,107)', width = 1.5))) %>%
add_trace(x = ~QuaterName, y = ~avg, type = 'scatter', mode = 'lines', name = 'Durchschnitt',
line = list(color = 'rgb(205, 12, 24)',width = 0.7),
hoverinfo = "text") %>%
layout(title = "Materialbedarf nach Quartalen",
bargap = 0.4,
xaxis = list(title = ""),
yaxis = list(title = ""),
legend = list(orientation = 'h'))
p
Dieser Plot zeigt die Verteilung der monatlichen Quantity Mengen nach Materialart an. Insbesondere A1662708602 hat starke Ausreißer nach oben.
theme_set(theme_ridges())
#Nach Date und Materialart gruppieren und summieren
df_p4 <- df_sendungen %>%
group_by(Date, Material) %>%
summarise(Quantity = sum(ShipmentQuantity))
#Die schönen Paletten haben nur 9 Farben -> 1 muss künstlich hinzugefügt werden
mycolors <- brewer.pal(n = 9, name = "Blues")
mycolors <- c(mycolors, "#001642")
ggplot(df_p4, aes(x = Quantity, y = Material)) +
geom_density_ridges(aes(fill = Material)) +
scale_fill_manual(values = mycolors)+
labs(title = 'Verteilung der monatlichen Quantity Mengen nach Materialart')
Dieser Plot zeigt die Struktur der zugrundeliegenden Supply Chain an. Dafür werden Techniken des Process Mining mit dem Bupar Package eingesetzt. Hierbei fungiert die ShipmentNumber als CaseID, der Hafen als Event und das Schiff als Resource. Hier ist zu erkennen wie die Vernetzung zwischen den Häfen ist.
df_pm <- data
#Die Aktivität wird definiert als Hafenbezeichnung
df_pm$Activity <- df_pm$Bezeichnung
#Nur diese Spalten selektieren
df_pm <- df_pm %>%
select(ShipmentNumber, Activity, VesselName, EventTime)
#Da die Daten auf der granularen Ebene der Bestellungen ist, müssen die Duplikate entfernt werden -> wir wollen nur die Sendungen verfolgen
df_pm <- unique(df_pm)
#Künstliche Activity Instance hinzufügen
df_pm <- df_pm %>%
mutate(activity_instance = 1:nrow(.))
#Künstlichen LifeCycle hinzufügen, da keine Informationen vorhanden sind für die Länge der Events
df_pm$LifeCycle <- "End"
#Duplikate entfernen
df_pm <- df_pm[!duplicated(df_pm[1:2]),]
#Event Log kreieren
event_logs <- df_pm %>%
eventlog(
case_id = "ShipmentNumber",
activity_id = "Activity",
activity_instance_id = "activity_instance",
lifecycle_id = "LifeCycle",
timestamp = "EventTime",
resource_id = "VesselName"
)
#Prozess Map plotten
event_logs %>%
filter_activity_frequency(percentage = .99) %>%
filter_trace_frequency(percentage = .99) %>%
process_map(type = frequency("absolute"))
event_logs %>%
filter_activity_frequency(percentage = .99) %>%
filter_trace_frequency(percentage = .99) %>%
process_map(processmapR::performance(mean, "days"))
event_logs %>%
precedence_matrix(type = "relative_consequent") %>%
plot
event_logs %>%
idle_time("resource", units = "days") %>%
plot()
Diese Weltkarte zeigt die Transportwege der Supply Chain. Die Größte der Punkte steht für die ankommende Menge an diesem Ort.
#Die Vorgänger Nachgänger Beziehungen extrahieren, um die Transportwege zu visualisieren
pred_matrix <- event_logs %>%
precedence_matrix(type = "relative_consequent")
#Diese Verbindungen haben keinen definierten Ort
df_map <- pred_matrix %>%
filter(consequent != "End") %>%
filter(antecedent != "Start")
#Daten joinen, um lat und long zu bekommen
df_map <- left_join(df_map, raw_event_location, by = c("antecedent" = "Bezeichnung"))
df_map <- df_map %>%
select(antecedent, consequent, rel_consequent, n, lat, long)
colnames(df_map) <- c("antecedent", "consequent", "rel_consequent", "n", "lat_antecedent", "long_antecedent")
df_map <- left_join(df_map, raw_event_location, by = c("consequent" = "Bezeichnung"))
df_map <- df_map %>%
select(antecedent, consequent, rel_consequent,n, lat_antecedent, long_antecedent, lat, long)
colnames(df_map) <- c("antecedent", "consequent", "rel_consequent","n", "lat_antecedent", "long_antecedent", "lat_consequent", "long_consequent")
#Zur Ansicht die Daten skalieren
df_map$nscale <- rescale(df_map$n, to=c(1,5))
#Size nach Anzahl der ankommenden Lieferungen und der ausgehenden Lieferungen -> Filtermöglichkeiten für später
ag <- aggregate(list(AnkommendeMenge=df_map$n), by=list(Bezeichnung=df_map$consequent), FUN=sum)
raw_event_location <- left_join(raw_event_location, ag, "Bezeichnung")
ag <- aggregate(list(AusgehendeMenge=df_map$n), by=list(Bezeichnung=df_map$antecedent), FUN=sum)
raw_event_location <- left_join(raw_event_location, ag, "Bezeichnung")
#NA mit 0 ersetzen
raw_event_location$AnkommendeMenge[is.na(raw_event_location$AnkommendeMenge)] <- 0
#Die Mengen skalieren, da sonst zu groß
raw_event_location$ScaledAnkommendeMenge <- rescale(raw_event_location$AnkommendeMenge, to=c(0.25,2))
raw_event_location$ScaledAusgehendeMenge <- rescale(raw_event_location$AusgehendeMenge, to=c(0.25,2))
#Fahrt benennen, um den hover Text anzupassen
df_map$FahrtName <- paste(df_map$antecedent,df_map$consequent, sep = " - ")
geo <- list(
scope = 'world',
projection = list(type = 'azimuthal equal area'),
showland = TRUE,
landcolor = toRGB("gray95"),
countrycolor = toRGB("gray80")
)
p <- plot_geo(locationmode = 'USA-states')
for(i in 1:nrow(raw_event_location)){
df_el <- as.data.frame(raw_event_location[i,])
p <- p %>%
add_markers(
data = df_el, x = ~long, y = ~lat, text = ~Bezeichnung,
size = ~ScaledAnkommendeMenge, hoverinfo = "text", alpha = 0.5
)
}
p <- p %>%
add_segments(
data = df_map,
x = ~long_antecedent, xend = ~long_consequent,
y = ~lat_antecedent, yend = ~lat_consequent, text = ~FahrtName,
alpha = 0.3, size = ~nscale, hoverinfo = "text"
) %>%
layout(
title = 'Supply Chain Routen - Mercedes Benz',
geo = geo, showlegend = FALSE
)
p
Nach reiflicher Überlegung wurden 7 Features ausgewählt, die einen potenziellen Einfluss auf die Ankuftszeit haben. Diese sind: -Material -EventName -EventLocation bzw. Bezeichnung -VesselName -Monat -Jahr -Wochentag
#Zuerst muss der Datensatz auf die ShipmentNumber reduziert werden -> Auf dieser Basis werden die Bestellungen rausgefiltert.
#Ansonsten hat man doppelte Daten drin
df <- data %>%
select(ShipmentNumber, Material, EventName, Bezeichnung, VesselName, monat, jahr, wochentag, ShipmentArrived, Time2Arrival)
df <- df[!duplicated(df[1:9]),]
#Herausfiltern von ShipmentArrived = TRUE. Macht keinen Sinn das forezucasten
df <- df %>%
filter(ShipmentArrived == FALSE)
#ShipmentArrived und ShipmentNumber droppen
df <- df %>%
select(-ShipmentArrived, -ShipmentNumber)
Mit One-Hot Encoding werden alle Categorial Variablen umgewandelt in numerische. Damit können die Machine Learning Verfahren umgehen.
colnames(df)
## [1] "Material" "EventName" "Bezeichnung" "VesselName"
## [5] "monat" "jahr" "wochentag" "Time2Arrival"
#Als Factor speichern
df$Material <- as.factor(df$Material)
df$EventName <- as.factor(df$EventName)
df$Bezeichnung <- as.factor(df$Bezeichnung)
df$VesselName <- as.factor(df$VesselName)
df$monat <- as.factor(df$monat)
df$jahr <- as.factor(df$jahr)
df$wochentag <- as.factor(df$wochentag)
#One-Hot Encoding für alle Spalten
df_tb <- setDT(df)
df_tb <- one_hot(df_tb, cols = "Material")
df_tb <- one_hot(df_tb, cols = "EventName")
df_tb <- one_hot(df_tb, cols = "Bezeichnung")
df_tb <- one_hot(df_tb, cols = "VesselName")
df_tb <- one_hot(df_tb, cols = "monat")
df_tb <- one_hot(df_tb, cols = "jahr")
df_tb <- one_hot(df_tb, cols = "wochentag")
#Wir haben jetzt 81 Features
df <- as.data.frame(df_tb)
colnames(df)
## [1] "Material_A1662707202"
## [2] "Material_A1662708602"
## [3] "Material_A1662708702"
## [4] "Material_A2132704801"
## [5] "Material_A2560107300"
## [6] "Material_A2572700000"
## [7] "Material_A2760101414"
## [8] "Material_A2760101614"
## [9] "Material_A2760106714"
## [10] "Material_A6540107406"
## [11] "EventName_Arrival Ship"
## [12] "EventName_Arrival Truck"
## [13] "EventName_Container Closed"
## [14] "EventName_Departure Ship"
## [15] "EventName_Departure Truck"
## [16] "EventName_Goods Receipt Dock"
## [17] "EventName_Loading Ship"
## [18] "EventName_Loading Truck"
## [19] "EventName_New Scheduling"
## [20] "EventName_Registration Yard"
## [21] "EventName_Unloading Ship"
## [22] "Bezeichnung_Bahnterminal Savannah"
## [23] "Bezeichnung_Consolidation Center Speyer"
## [24] "Bezeichnung_Container Dummy"
## [25] "Bezeichnung_Hafen Aguascalientes"
## [26] "Bezeichnung_Hafen Altamira"
## [27] "Bezeichnung_Hafen Antwerpen"
## [28] "Bezeichnung_Hafen Bremerhaven"
## [29] "Bezeichnung_Hafen Charleston"
## [30] "Bezeichnung_Hafen Hamburg"
## [31] "Bezeichnung_Hafen Rotterdam"
## [32] "Bezeichnung_Werk Tuscaloosa"
## [33] "VesselName_BREVIK BRIDGE"
## [34] "VesselName_CHARLESTON EXPRESS"
## [35] "VesselName_EASLINE DALIAN"
## [36] "VesselName_MAINE TRADER"
## [37] "VesselName_MOL EMISSARY"
## [38] "VesselName_MOL EMPIRE"
## [39] "VesselName_MOL EXPERIENCE"
## [40] "VesselName_MONACO MAERSK"
## [41] "VesselName_MUNICH MAERSK"
## [42] "VesselName_NYK DAEDALUS"
## [43] "VesselName_NYK DENEB"
## [44] "VesselName_NYK METEOR"
## [45] "VesselName_NYK NEBULA"
## [46] "VesselName_NYK REMUS"
## [47] "VesselName_NYK RIGEL"
## [48] "VesselName_NYK ROMULUS"
## [49] "VesselName_NYK RUMINA"
## [50] "VesselName_PHILADELPHIA EXPRESS"
## [51] "VesselName_ROTTERDAM EXPRESS"
## [52] "VesselName_SHANGHAI TRADER"
## [53] "VesselName_ST LOUIS EXPRESS"
## [54] "VesselName_TRF PARTICI"
## [55] "VesselName_VECCHIO BRIDGE"
## [56] "VesselName_VENICE BRIDGE"
## [57] "VesselName_WASHINGTON EXPRESS"
## [58] "VesselName_YM UNITY"
## [59] "VesselName_YORKTOWN EXPRESS"
## [60] "monat_1"
## [61] "monat_2"
## [62] "monat_3"
## [63] "monat_4"
## [64] "monat_5"
## [65] "monat_6"
## [66] "monat_7"
## [67] "monat_8"
## [68] "monat_9"
## [69] "monat_10"
## [70] "monat_11"
## [71] "monat_12"
## [72] "jahr_2018"
## [73] "jahr_2019"
## [74] "wochentag_Dienstag"
## [75] "wochentag_Donnerstag"
## [76] "wochentag_Freitag"
## [77] "wochentag_Mittwoch"
## [78] "wochentag_Montag"
## [79] "wochentag_Samstag"
## [80] "wochentag_Sonntag"
## [81] "Time2Arrival"
remove(df_tb)
Um das Modell anschließend evaluieren zu können, wird der Datensatz unterteilt in einen Trainingsset und Testset (80/20)
#Definition der Sample size
smp_size <- floor(0.8 * nrow(df))
#Set Seed
set.seed(100)
train_ind <- sample(seq_len(nrow(df)), size = smp_size)
train <- df[train_ind, ]
test <- df[-train_ind, ]
Da relativ wenige Daten vorhanden sind, werden voraussichtlich neuronale Netze nicht die besten Ergebnisse liefern. Sie werden hier trotzdem mit aufgeführt, um für die Skalierbarkeit bereits das Modell zu haben.
#Als Datatable speichern
train_x <- setDT(train)
test_x <- setDT(test)
#Targetvariablen speichern
labels_train <- train_x$Time2Arrival
labels_test <- test_x$Time2Arrival
#Targetvariablen entfernen, sodass nur Features da sind
train_x <- subset(train_x, select = - Time2Arrival)
test_x <- subset(test_x, select = - Time2Arrival)
#Als Matrix speichern
train_x <- data.matrix(train_x, rownames.force = NA)
test_x <- data.matrix(test_x, rownames.force = NA)
#Als XGBoost MAtrix speichern
dtrain <- xgb.DMatrix(data = train_x,label = labels_train)
dtest <- xgb.DMatrix(data = test_x,label=labels_test)
#Für die Evaluation während des Trainings
watchlist <- list(train=dtrain, test=dtest)
set.seed(101)
st <- xgb.train(data=dtrain, booster = "gbtree", max.depth=7, lambda = 2, alpha = 2, colsample_bytree = 0.5, min_child_weight = 0.2, early_stopping_rounds= 200, nthread = 2, gamma = 0, eta= 0.01, nrounds=10000, watchlist=watchlist, eval.metric = "rmse", objective = "reg:linear")
#Auf die Testdaten predicten
xgbpred <- predict (st,dtest)
#Für späteres Stacking
xgbpred_train <- predict(st, dtrain)
#xgb.save(st, "xgb.model")
Das vorher erstellte Modell wird nun anhand der Testdaten validiert.
#Validation Kennzahlen
rmse(labels_test, round(xgbpred,0))
## [1] 17.62393
MAE(labels_test, round(xgbpred,0))
## [1] 9.39945
#Feautre Importance
mat <- xgb.importance(feature_names = colnames(train_x),model = st)
xgb.plot.importance (importance_matrix = mat[1:30])
#Ausschnitt des Tree's
model <- xgb.dump(st, with.stats = T)
model[1:30]
## [1] "booster[0]"
## [2] "0:[f19<0.5] yes=1,no=2,missing=1,gain=1348963,cover=20366"
## [3] "1:[f28<0.5] yes=3,no=4,missing=3,gain=848633,cover=17549"
## [4] "3:[f71<0.5] yes=7,no=8,missing=7,gain=290917,cover=8785"
## [5] "7:[f21<0.5] yes=15,no=16,missing=15,gain=33451.75,cover=1980"
## [6] "15:[f5<0.5] yes=31,no=32,missing=31,gain=6525.875,cover=1784"
## [7] "31:[f3<0.5] yes=57,no=58,missing=57,gain=2595.375,cover=1537"
## [8] "57:[f33<0.5] yes=93,no=94,missing=93,gain=1005.375,cover=1321"
## [9] "93:leaf=0.280295104,cover=1184"
## [10] "94:leaf=0.327302158,cover=137"
## [11] "58:leaf=0.237935767,cover=216"
## [12] "32:leaf=0.217811242,cover=247"
## [13] "16:[f45<0.5] yes=33,no=34,missing=33,gain=1063.27344,cover=196"
## [14] "33:[f5<0.5] yes=59,no=60,missing=59,gain=2070.19922,cover=183"
## [15] "59:[f46<0.5] yes=95,no=96,missing=95,gain=173.857422,cover=177"
## [16] "95:leaf=0.133896098,cover=152"
## [17] "96:leaf=0.0850000009,cover=25"
## [18] "60:leaf=0.308750004,cover=6"
## [19] "34:leaf=0.0350000001,cover=13"
## [20] "8:[f21<0.5] yes=17,no=18,missing=17,gain=158137,cover=6805"
## [21] "17:[f30<0.5] yes=35,no=36,missing=35,gain=124457,cover=6551"
## [22] "35:[f5<0.5] yes=61,no=62,missing=61,gain=31794,cover=6422"
## [23] "61:[f3<0.5] yes=97,no=98,missing=97,gain=31640,cover=5675"
## [24] "97:leaf=0.41574043,cover=5049"
## [25] "98:leaf=0.336449057,cover=626"
## [26] "62:[f29<0.5] yes=99,no=100,missing=99,gain=5413.625,cover=747"
## [27] "99:leaf=0.344053268,cover=674"
## [28] "100:leaf=0.236733332,cover=73"
## [29] "36:[f60<0.5] yes=63,no=64,missing=63,gain=36895.6875,cover=129"
## [30] "63:leaf=0.599607825,cover=100"
train_labels <- labels_train
train_data <- as.data.frame(train_x)
test_labels <- labels_test
test_data <- as.data.frame(test_x)
#Dataframe als numeric speichern
train_data <- as.data.frame(lapply(train_data, as.numeric))
test_data <- as.data.frame(lapply(test_data, as.numeric))
#Check ob wirklich alle NA's raus sind -> Neuronale Netze können damit nicht umgehen
colnames(train_data)[colSums(is.na(train_data)) > 0]
## character(0)
colnames(test_data)[colSums(is.na(test_data)) > 0]
## character(0)
#Als Matrix speichern
train_data <- as.matrix(train_data)
test_data <- as.matrix(test_data)
#Neuronales Netz mit einem Input Layer, ein Hidden Layer und ein Output Layer
build_model <- function() {
model <- keras_model_sequential() %>%
layer_dense(units = 64, activation = "relu", input_shape = dim(train_data)[2]) %>%
layer_dense(units = 32, activation = "relu") %>%
layer_dense(units = 1)
model %>% compile(
loss = 'mean_squared_error',
optimizer = optimizer_rmsprop(),
metrics = list("mean_absolute_error")
)
model
}
model <- build_model()
model %>% summary()
## ___________________________________________________________________________
## Layer (type) Output Shape Param #
## ===========================================================================
## dense (Dense) (None, 64) 5184
## ___________________________________________________________________________
## dense_1 (Dense) (None, 32) 2080
## ___________________________________________________________________________
## dense_2 (Dense) (None, 1) 33
## ===========================================================================
## Total params: 7,297
## Trainable params: 7,297
## Non-trainable params: 0
## ___________________________________________________________________________
#Funktion damit das Modell die Scores während des Trainings ausgibt
print_dot_callback <- callback_lambda(
on_epoch_end = function(epoch, logs) {
if (epoch %% 80 == 0) cat("\n")
cat(".")
}
)
#Rundenanzahl
epochs <- 500
#Early Stopping wenn sich der Score nicht verbessert
early_stop <- callback_early_stopping(monitor = "val_loss", patience = 20)
model <- build_model()
history <- model %>% fit(
train_data,
train_labels,
epochs = epochs,
validation_split = 0.2,
verbose = 0,
batch_size = 32,
callbacks = list(early_stop, print_dot_callback)
)
##
## ................................................................................
## ........................
test_predictions <- model %>% predict(test_data)
train_predictions_nn <- model %>% predict(train_data)
Das vorher erstellte Modell wird nun anhand der Testdaten validiert.
#Validation Kennzahlen
rmse(test_labels, test_predictions[ , 1])
## [1] 17.68879
MAE(test_labels, test_predictions[ , 1])
## [1] 9.293606
plot(history, metrics = "mean_absolute_error", smooth = FALSE)
linearMod <- lm(Time2Arrival ~. , data=train)
predictions_lm = predict.lm(linearMod, test)
prediction_lm_train <- predict.lm(linearMod, train)
rmse(test_labels, predictions_lm)
## [1] 20.18799
MAE(test_labels, predictions_lm)
## [1] 11.48298
summary(linearMod)
##
## Call:
## lm(formula = Time2Arrival ~ ., data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -75.522 -9.842 -3.287 3.770 228.157
##
## Coefficients: (11 not defined because of singularities)
## Estimate Std. Error t value
## (Intercept) -12.53703 3.52240 -3.559
## Material_A1662707202 2.14110 0.65749 3.256
## Material_A1662708602 1.57064 0.73294 2.143
## Material_A1662708702 1.95924 0.65278 3.001
## Material_A2132704801 -0.25479 0.68912 -0.370
## Material_A2560107300 9.67005 0.74869 12.916
## Material_A2572700000 -1.52356 0.67103 -2.270
## Material_A2760101414 3.61161 0.76921 4.695
## Material_A2760101614 7.95987 0.72843 10.927
## Material_A2760106714 2.96604 0.76996 3.852
## Material_A6540107406 NA NA NA
## `EventName_Arrival Ship` -3.50613 0.54633 -6.418
## `EventName_Arrival Truck` 46.06577 2.75694 16.709
## `EventName_Container Closed` 31.31717 4.53156 6.911
## `EventName_Departure Ship` 39.70288 2.80878 14.135
## `EventName_Departure Truck` -1.93770 0.52535 -3.688
## `EventName_Goods Receipt Dock` NA NA NA
## `EventName_Loading Ship` 40.15621 2.80905 14.295
## `EventName_Loading Truck` 0.69660 4.42213 0.158
## `EventName_New Scheduling` -0.78814 0.77825 -1.013
## `EventName_Registration Yard` -10.28634 3.30757 -3.110
## `EventName_Unloading Ship` NA NA NA
## `Bezeichnung_Bahnterminal Savannah` -4.55551 3.50357 -1.300
## `Bezeichnung_Consolidation Center Speyer` NA NA NA
## `Bezeichnung_Container Dummy` NA NA NA
## `Bezeichnung_Hafen Aguascalientes` 56.78961 19.97288 2.843
## `Bezeichnung_Hafen Altamira` 105.34659 14.44532 7.293
## `Bezeichnung_Hafen Antwerpen` -35.77509 2.25407 -15.871
## `Bezeichnung_Hafen Bremerhaven` -19.01189 1.84825 -10.286
## `Bezeichnung_Hafen Charleston` 7.74333 3.26426 2.372
## `Bezeichnung_Hafen Hamburg` -23.71123 1.91485 -12.383
## `Bezeichnung_Hafen Rotterdam` NA NA NA
## `Bezeichnung_Werk Tuscaloosa` NA NA NA
## `VesselName_BREVIK BRIDGE` 6.77549 2.62423 2.582
## `VesselName_CHARLESTON EXPRESS` 2.11874 0.64729 3.273
## `VesselName_EASLINE DALIAN` 29.09705 4.04302 7.197
## `VesselName_MAINE TRADER` -3.15182 1.04452 -3.017
## `VesselName_MOL EMISSARY` 2.45330 0.66488 3.690
## `VesselName_MOL EMPIRE` 0.14822 0.74267 0.200
## `VesselName_MOL EXPERIENCE` 9.34574 2.15481 4.337
## `VesselName_MONACO MAERSK` -1.06983 1.83745 -0.582
## `VesselName_MUNICH MAERSK` -1.37810 3.09580 -0.445
## `VesselName_NYK DAEDALUS` 9.99980 3.11126 3.214
## `VesselName_NYK DENEB` 13.67662 4.72397 2.895
## `VesselName_NYK METEOR` 22.23586 4.22302 5.265
## `VesselName_NYK NEBULA` 14.94052 4.26605 3.502
## `VesselName_NYK REMUS` 7.65935 4.13603 1.852
## `VesselName_NYK RIGEL` 5.15349 2.79786 1.842
## `VesselName_NYK ROMULUS` 10.59590 3.37852 3.136
## `VesselName_NYK RUMINA` 10.43309 1.44124 7.239
## `VesselName_PHILADELPHIA EXPRESS` 0.24397 0.61269 0.398
## `VesselName_ROTTERDAM EXPRESS` -5.99036 1.37727 -4.349
## `VesselName_SHANGHAI TRADER` -1.25716 0.68483 -1.836
## `VesselName_ST LOUIS EXPRESS` 0.52678 0.66818 0.788
## `VesselName_TRF PARTICI` -0.09635 1.06666 -0.090
## `VesselName_VECCHIO BRIDGE` -2.23403 0.74827 -2.986
## `VesselName_VENICE BRIDGE` 10.57486 2.34041 4.518
## `VesselName_WASHINGTON EXPRESS` -0.52700 0.61570 -0.856
## `VesselName_YM UNITY` 0.27787 5.80456 0.048
## `VesselName_YORKTOWN EXPRESS` NA NA NA
## monat_1 24.72397 1.14191 21.651
## monat_2 23.26608 1.08257 21.492
## monat_3 17.17387 0.95939 17.901
## monat_4 14.94688 0.90336 16.546
## monat_5 9.07486 1.04600 8.676
## monat_6 0.23413 0.88662 0.264
## monat_7 -0.88576 0.90192 -0.982
## monat_8 3.35112 0.86115 3.891
## monat_9 1.35695 0.85489 1.587
## monat_10 4.04854 0.86612 4.674
## monat_11 6.21375 0.96150 6.463
## monat_12 NA NA NA
## jahr_2018 26.26187 0.81802 32.104
## jahr_2019 NA NA NA
## wochentag_Dienstag -3.52513 0.56717 -6.215
## wochentag_Donnerstag -2.60314 0.62327 -4.177
## wochentag_Freitag -4.95353 0.91097 -5.438
## wochentag_Mittwoch -3.48288 0.59836 -5.821
## wochentag_Montag -2.70880 0.52936 -5.117
## wochentag_Samstag -5.25146 0.90564 -5.799
## wochentag_Sonntag NA NA NA
## Pr(>|t|)
## (Intercept) 0.000373 ***
## Material_A1662707202 0.001130 **
## Material_A1662708602 0.032131 *
## Material_A1662708702 0.002691 **
## Material_A2132704801 0.711589
## Material_A2560107300 < 2e-16 ***
## Material_A2572700000 0.023188 *
## Material_A2760101414 2.68e-06 ***
## Material_A2760101614 < 2e-16 ***
## Material_A2760106714 0.000117 ***
## Material_A6540107406 NA
## `EventName_Arrival Ship` 1.41e-10 ***
## `EventName_Arrival Truck` < 2e-16 ***
## `EventName_Container Closed` 4.96e-12 ***
## `EventName_Departure Ship` < 2e-16 ***
## `EventName_Departure Truck` 0.000226 ***
## `EventName_Goods Receipt Dock` NA
## `EventName_Loading Ship` < 2e-16 ***
## `EventName_Loading Truck` 0.874832
## `EventName_New Scheduling` 0.311210
## `EventName_Registration Yard` 0.001874 **
## `EventName_Unloading Ship` NA
## `Bezeichnung_Bahnterminal Savannah` 0.193531
## `Bezeichnung_Consolidation Center Speyer` NA
## `Bezeichnung_Container Dummy` NA
## `Bezeichnung_Hafen Aguascalientes` 0.004469 **
## `Bezeichnung_Hafen Altamira` 3.15e-13 ***
## `Bezeichnung_Hafen Antwerpen` < 2e-16 ***
## `Bezeichnung_Hafen Bremerhaven` < 2e-16 ***
## `Bezeichnung_Hafen Charleston` 0.017694 *
## `Bezeichnung_Hafen Hamburg` < 2e-16 ***
## `Bezeichnung_Hafen Rotterdam` NA
## `Bezeichnung_Werk Tuscaloosa` NA
## `VesselName_BREVIK BRIDGE` 0.009833 **
## `VesselName_CHARLESTON EXPRESS` 0.001065 **
## `VesselName_EASLINE DALIAN` 6.38e-13 ***
## `VesselName_MAINE TRADER` 0.002552 **
## `VesselName_MOL EMISSARY` 0.000225 ***
## `VesselName_MOL EMPIRE` 0.841808
## `VesselName_MOL EXPERIENCE` 1.45e-05 ***
## `VesselName_MONACO MAERSK` 0.560412
## `VesselName_MUNICH MAERSK` 0.656215
## `VesselName_NYK DAEDALUS` 0.001311 **
## `VesselName_NYK DENEB` 0.003794 **
## `VesselName_NYK METEOR` 1.41e-07 ***
## `VesselName_NYK NEBULA` 0.000462 ***
## `VesselName_NYK REMUS` 0.064060 .
## `VesselName_NYK RIGEL` 0.065498 .
## `VesselName_NYK ROMULUS` 0.001714 **
## `VesselName_NYK RUMINA` 4.68e-13 ***
## `VesselName_PHILADELPHIA EXPRESS` 0.690492
## `VesselName_ROTTERDAM EXPRESS` 1.37e-05 ***
## `VesselName_SHANGHAI TRADER` 0.066414 .
## `VesselName_ST LOUIS EXPRESS` 0.430480
## `VesselName_TRF PARTICI` 0.928023
## `VesselName_VECCHIO BRIDGE` 0.002834 **
## `VesselName_VENICE BRIDGE` 6.27e-06 ***
## `VesselName_WASHINGTON EXPRESS` 0.392039
## `VesselName_YM UNITY` 0.961820
## `VesselName_YORKTOWN EXPRESS` NA
## monat_1 < 2e-16 ***
## monat_2 < 2e-16 ***
## monat_3 < 2e-16 ***
## monat_4 < 2e-16 ***
## monat_5 < 2e-16 ***
## monat_6 0.791725
## monat_7 0.326072
## monat_8 1.00e-04 ***
## monat_9 0.112466
## monat_10 2.97e-06 ***
## monat_11 1.05e-10 ***
## monat_12 NA
## jahr_2018 < 2e-16 ***
## jahr_2019 NA
## wochentag_Dienstag 5.22e-10 ***
## wochentag_Donnerstag 2.97e-05 ***
## wochentag_Freitag 5.46e-08 ***
## wochentag_Mittwoch 5.95e-09 ***
## wochentag_Montag 3.13e-07 ***
## wochentag_Samstag 6.78e-09 ***
## wochentag_Sonntag NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.86 on 20296 degrees of freedom
## Multiple R-squared: 0.3584, Adjusted R-squared: 0.3563
## F-statistic: 164.3 on 69 and 20296 DF, p-value: < 2.2e-16
#Dient mehr als Pilot
#Die Ergebnisse der vorherigen Modelle werden gestacked. Darauf wird anschließend eine Support Vector Regression angewendet.
stacked_train <- as.data.frame(xgbpred_train)
stacked_train$nn <- train_predictions_nn[ , 1]
stacked_train$lr <- prediction_lm_train
colnames(stacked_train) <- c("xgb","nn","lr")
stacked_train$target <- train$Time2Arrival
stacked_test <- as.data.frame(xgbpred)
stacked_test$nn <- test_predictions[ , 1]
stacked_test$lr <- predictions_lm
colnames(stacked_test) <- c("xgb","nn","lr")
#Dieser Plot stellt die Predictions der beiden Modelle gegenüber.
plot_ly(data = stacked_train, x = ~xgb, y = ~nn) %>%
layout(title = "XGB vs. NN")
set.seed(1005)
#SVM
svm_model <- svm(target ~. , stacked_train)
test_svm <- predict(svm_model, stacked_test)
rmse(test_labels, test_svm)
## [1] 17.3852
MAE(test_labels, test_svm)
## [1] 8.791122
#Ursprungsdatensatz mit normalen Features
df_check <- data %>%
select(ShipmentNumber, Material, EventName, Bezeichnung, VesselName, monat, jahr, wochentag, ShipmentArrived, Time2Arrival)
df_check <- df_check [!duplicated(df_check [1:9]),]
df_check <- df_check %>%
filter(ShipmentArrived == FALSE)
#Gleichen Datensatz wiederherstellen
test_check <- df_check[-train_ind, ]
#Prediction hinzufügen
test_check$prediction <- round(test_svm,0)
#Absoluter Error
test_check$error <- abs(test_check$Time2Arrival - test_check$prediction)
#Nach Häfen gruppieren und durchschnittlichen Error berechnen
df_check <- test_check %>%
group_by(Bezeichnung) %>%
summarise(AVG = mean(error))
df_check$AVG <- round(df_check$AVG,0)
plot_ly(df_check, x = ~Bezeichnung, y = ~AVG, type = 'bar', text = ~AVG, textposition = 'auto',
marker = list(color = 'rgb(158,202,225)',
line = list(color = 'rgb(8,48,107)', width = 1.5))) %>%
layout(title = "Mean Error nach Häfen",
xaxis = list(title = ""),
yaxis = list(title = ""))
Die bisherigen Ergebnisse sind sehr vielversprechend. Das weitere Vorgehen wird in Abstimmung mit dem Management und den Verantwortlichen in den Fachbereichen diskutiert.